version 10.0
log using "C:\fsu2010\formateur\replication\changepm.log", replace
#delimit ;
set mem 30m;
/*clear matrix;*/
/*clear mata; */

set more off;

*     ***************************************************************** *;
*     ***************************************************************** *;
*       File-Name:  changepm.do                                         *;
*       Date:       November 9, 2010                                    *;
*       Author:     GG/MG                                               *;
*       Purpose:    Produce Table 5b in AJPS paper                      *;
*       Input File:     europe.dta                                      *;
*       Output File:    changepm.log                                    *;
*       Data Output:    none                                            *;
*       Previous file:  europe.dta                                      *;
*       Machine:        Matt's desktop                                  *;
*     ****************************************************************  *;
*     ****************************************************************  *;

use "C:\fsu2010\formateur\replication\europe.dta", clear;


*     ****************************************************************  *;
*       Model 6: Western Europe                                         *;
*     ****************************************************************  *;

mixlogit choice largestparty_east party_seatshare_east medianparty1_east presidentparty_east previousPM_east previousPM_conflict_east previousPM_gain1_east previous_cabinet_gain1_east,  
    rand (largestparty party_seatshare medianparty1 presidentparty previousPM president_investiture president_direct previousPM_conflict previousPM_gain1 previous_cabinet_gain1) 
    ln(0) nrep(125) group(cabinetcode);
    
*     ****************************************************************  *;
*       Now produce Figure 1                                            *;
*     ****************************************************************  *;

matrix coeffs = e(b);
matrix covmat = e(V);

keep if cabinetcode==1222;

set obs 10000;

generate pvda = 0;
replace pvda = 1 if partyID==7;
generate cda = 0;
replace cda = 1 if partyID==13;
generate vvd = 0;
replace vvd = 1 if partyID==22;
generate d66 = 0;
replace d66 = 1 if partyID==10;

generate ppvda1 = .;
generate pcda1 = .;
generate pvvd1 = .;
generate pd661 = .;
generate ppvda2 = .;
generate pcda2 = .;
generate pvvd2 = .;
generate pd662 = .;
generate ppvda3 = .;
generate pcda3 = .;
generate pvvd3 = .;
generate pd663 = .;
generate ppvda4 = .;
generate pcda4 = .;
generate pvvd4 = .;
generate pd664 = .;

drawnorm b1-b28, means(coeffs) cov(covmat) double; 

nois _dots 0, title(largest party sim) reps(10000);
forvalues j = 1/10000 {;

  mkmat b1-b28 if _n==`j', matrix(bdraw);
  
mixlpred2 pprobs; 

  qui generate pt1 = pprobs * pvda;
  qui generate ct1 = pprobs * cda;
  qui generate vt1 = pprobs * vvd;
  qui generate dt1 = pprobs * d66;

  qui replace previousPM = 1 if pvda==1;
  qui replace previousPM_gain1 = -13.33333 if pvda==1;
  qui replace previous_cabinet_gain1 = -7.333344 if pvda==1;
  qui replace previousPM = 0 if cda==1;
  qui replace previousPM_gain1 = 0 if cda==1;
  qui replace previous_cabinet_gain1 = 0 if cda==1;

mixlpred2 pprobs2;

  qui generate pt2 = pprobs2 * pvda;
  qui generate ct2 = pprobs2 * cda;
  qui generate vt2 = pprobs2 * vvd;
  qui generate dt2 = pprobs2 * d66;

  qui replace previousPM = 1 if vvd==1;
  qui replace previousPM_gain1 = -13.33333 if vvd==1;
  qui replace previous_cabinet_gain1 = -7.333344 if vvd==1;
  qui replace previousPM = 0 if pvda==1;
  qui replace previousPM_gain1 = 0 if pvda==1;
  qui replace previous_cabinet_gain1 = 0 if pvda==1;

mixlpred2 pprobs3;

  qui generate pt3 = pprobs3 * pvda;
  qui generate ct3 = pprobs3 * cda;
  qui generate vt3 = pprobs3 * vvd;
  qui generate dt3 = pprobs3 * d66;

  qui replace previousPM = 1 if d66==1;
  qui replace previousPM_gain1 = -13.33333 if d66==1;
  qui replace previous_cabinet_gain1 = -7.333344 if d66==1;
  qui replace previousPM = 0 if vvd==1;
  qui replace previousPM_gain1 = 0 if vvd==1;
  qui replace previous_cabinet_gain1 = 0 if vvd==1;

mixlpred2 pprobs4;

  qui generate pt4 = pprobs4 * pvda;
  qui generate ct4 = pprobs4 * cda;
  qui generate vt4 = pprobs4 * vvd;
  qui generate dt4 = pprobs4 * d66;

  qui replace previousPM = 1 if cda==1;
  qui replace previousPM_gain1 = -13.33333 if cda==1;
  qui replace previous_cabinet_gain1 = -7.333344 if cda==1;
  qui replace previousPM = 0 if d66==1;
  qui replace previousPM_gain1 = 0 if d66==1;
  qui replace previous_cabinet_gain1 = 0 if d66==1;


  egen pt1s = sum(pt1);
  egen ct1s = sum(ct1);
  egen vt1s = sum(vt1);
  egen dt1s = sum(dt1);
  egen pt2s = sum(pt2);
  egen ct2s = sum(ct2);
  egen vt2s = sum(vt2);
  egen dt2s = sum(dt2);
  egen pt3s = sum(pt3);
  egen ct3s = sum(ct3);
  egen vt3s = sum(vt3);
  egen dt3s = sum(dt3);
  egen pt4s = sum(pt4);
  egen ct4s = sum(ct4);
  egen vt4s = sum(vt4);
  egen dt4s = sum(dt4);

  qui replace ppvda1 = pt1s if _n==`j';
  qui replace pcda1 = ct1s if _n==`j';
  qui replace pvvd1 = vt1s if _n==`j';
  qui replace pd661 = dt1s if _n==`j';
  qui replace ppvda2 = pt2s if _n==`j';
  qui replace pcda2 = ct2s if _n==`j';
  qui replace pvvd2 = vt2s if _n==`j';
  qui replace pd662 = dt2s if _n==`j';
  qui replace ppvda3 = pt3s if _n==`j';
  qui replace pcda3 = ct3s if _n==`j';
  qui replace pvvd3 = vt3s if _n==`j';
  qui replace pd663 = dt3s if _n==`j';
  qui replace ppvda4 = pt4s if _n==`j';
  qui replace pcda4 = ct4s if _n==`j';
  qui replace pvvd4 = vt4s if _n==`j';
  qui replace pd664 = dt4s if _n==`j';

drop pt1 ct1 vt1 dt1 pt2 ct2 vt2 dt2 pt3 ct3 vt3 dt3 pt4 ct4 vt4 dt4 
pt1s ct1s vt1s dt1s pt2s ct2s vt2s dt2s pt3s ct3s vt3s dt3s pt4s ct4s vt4s dt4s 
pprobs pprobs2 pprobs3 pprobs4; 

nois _dots `j' 0;
  
};

generate pvdad2 = ppvda2 - ppvda1;
generate pvdad3 = ppvda3 - ppvda1;
generate pvdad4 = ppvda4 - ppvda1;
generate cdad2 = pcda2 - pcda1;
generate cdad3 = pcda3 - pcda1;
generate cdad4 = pcda4 - pcda1;
generate vvdd2 = pvvd2 - pvvd1;
generate vvdd3 = pvvd3 - pvvd1;
generate vvdd4 = pvvd4 - pvvd1;
generate d66d2 = pd662 - pd661;
generate d66d3 = pd663 - pd661;
generate d66d4 = pd664 - pd661;

summarize 
ppvda1 pcda1 pvvd1 pd661
ppvda2 pcda2 pvvd2 pd662
ppvda3 pcda3 pvvd3 pd663
ppvda4 pcda4 pvvd4 pd664
pvdad2 cdad2 vvdd2 d66d2
pvdad3 cdad3 vvdd3 d66d3
pvdad4 cdad4 vvdd4 d66d4;
